//Kernel kmeans.
//Kn: Normalised Kernel Matrix


inline
k_kmeans::k_kmeans(const mat &in_Kn, const int in_Ncent)
: K(in_Kn), N_cent(in_Ncent), N_data(K.n_cols)
{
}

inline
void
k_kmeans::calc(const uword max_iter)
{
  
  uword iter = 0;
  umat comparing;
  init_partitions();
  
  //while ( err > 0.005 && iter < max_iter)
  do 
  {  
    cout << "Iter= "  << iter;
    //update_kkmeans();
    update_kkmeans2();//-- Nov15,2013 
    comparing = (partitions != new_partitions);
    //cout << "See in Matlab" << endl;
    //partitions.save("old_partitions.mat", raw_ascii);
    //new_partitions.save("new_partitions.mat", raw_ascii);
    //getchar();
    //partitions.print("partitions");
    //comparing.print("comparing: ");
    if (iter < max_iter)
    {
      err = sum(sum(comparing));
      //err = 1;
      
    }
    else
    {
      err = 0;
      distortion = 0; 
    }
    //counts.t().print("counts");
    //err = 1;
    cout << "- Err= " << err << " & distortion= " << distortion << endl;
    partitions = new_partitions;
    counts = new_counts;
    iter++;
    //cout << "(err !=0 )  || (iter < max_iter) = " << ((err !=0 )  || (iter < max_iter)) << endl;
  } while ( (err !=0 )  );
  
  //cout << "Saving Lables: " << endl;
  //save_labels("labels_grass.dat");
  //partitions.save("partitions.mat", raw_ascii);
  //counts.t().print("Final counts");
  //Points_cent.save("min_distances.mat", raw_ascii);
}


inline
void
k_kmeans::init_partitions()
{
  cout << "Ndata: " << N_data << endl;  
  cout << "Ncent: " << N_cent << endl;
  
  //aux_functions auxfx;
  vec tmp = randu<vec>(N_data ) * (N_cent);
  //tmp.t().print("tmp:");
  uvec tmp2 = conv_to<uvec>::from( tmp );
  //tmp2.t().print("tmp2:");
  
  partitions.zeros(N_data,N_cent);
  new_partitions.zeros(N_data,N_cent);
  
  counts.zeros(N_cent);
  
  for (uword c = 0 ; c < N_cent; ++c)
  {
    uvec q1 = find (tmp2 == c);
    counts(c) = q1.n_elem;
    //q1.print("q1 bf: ");
    q1.resize(N_data);
    //q1.print("q1 after: ");
    partitions.col( c ) = q1;
  }
  //partitions.print("initial_indices:");
  partitions.save("ini_partitions.mat", raw_ascii);
  //uvec counts_sort  = sort(counts);
  counts.t().print("Ini counts:");
  //getchar();
  
}



inline
void
k_kmeans::update_kkmeans()
{
  vec dist(N_cent);
  double dist_cent;
  uword best_ci;
  distortion = 0;
  
  Points_cent.zeros(N_data,N_cent);// Nov14,2013 Same as partitions but with distances
  
  new_counts.zeros(N_cent);
  new_partitions.zeros(N_data,N_cent);
  
  for (uword ni = 0; ni < N_data; ++ni)
  {
    //cout << "Doing for ni = " << ni<< endl;
    for (uword c = 0; c < N_cent; ++c)
    { 
      cout << "Doing for cluster = " << c << "...";
      cout << "counts(c)= " << counts(c) << endl;
      getchar();
      double sum1 = 0;
      for (uword i = 0; i < counts(c); ++i)
      {
	//cout << "K(ni,partitions(i,c))= " << K(ni,partitions(i,c)) << endl;
	sum1 += K(ni,partitions(i,c));
      }
      
      //cout << "sum1= " << sum1 << endl;
      //getchar();
      double sum2 = 0;
      
      for (uword i = 0; i < counts(c); ++i)
      {
	for (uword j = 0; j < counts(c); ++j)
	{
	  sum2 += K(partitions(i,c), partitions(j,c));
	}
      }
      dist_cent = K(ni,ni) -2*sum1/counts(c) + sum2 / pow(counts(c),2 );
      dist(c) = dist_cent;
      //cout << "dist_cent: " << dist_cent << endl;
      
      
    }
    
    
    
    //dist.t().print("dist: ");
    double mind = dist.min(best_ci);
    distortion+=mind;
    //cout << "best_ci= " << best_ci << endl;  
    //getchar();
    new_partitions(new_counts(best_ci),best_ci) = ni;
    Points_cent (new_counts(best_ci),best_ci) = mind;
    new_counts(best_ci) += 1;
    
    //new_partitions.print("New partition: ");
    //new_counts.t().print("New counts: ");
    
  }    
  //new_counts.t().print("new_counts: ");
  //getchar();
  //cout << "press a key " << endl;
  //getchar();
}

//This function is only useful to plot
//the labels in Matlab.
inline
void
k_kmeans::save_labels(std::string name)
{
  uvec labels;
  labels.zeros(N_data);
  // uvec q1 = find (tmp2 == c);
  
  for (uword ci = 0; ci < N_cent; ++ci)  
  {
    uvec v = partitions.col(ci);
    //v.print("v: ");
    v.resize(counts(ci));
    //v.print("v: ");
    //labels(v).fill(ci);
    for (uword n = 0; n< counts(ci); ++n)
    {
      labels(v(n)) = ci;
    }
    
    //labels.print("labels: ");
    //getchar();*/
  }
  //labels.t().print("labels: ");
  labels.save(name, raw_ascii);
}


inline
void
k_kmeans::print_labels(const field<std::string>  labels)
{
  //partitions.print("Final Partitions");
  field<std::string> lab_partitions( N_cent, N_data  );
  // uvec q1 = find (tmp2 == c);
  
  for (uword ci = 0; ci < N_cent; ++ci)  
  {
    uvec v = partitions.col(ci);
    //v.print("v: ");
    v.resize(counts(ci));
    //v.print("v: ");
    //labels(v).fill(ci);
    cout << "Cluster " << ci + 1 << " has  " << counts(ci) << " points" << endl;
    for (uword n = 0; n< counts(ci); ++n)
    {
      //labels(v(n),0) = lab(ci,0);
      cout <<  labels( partitions( n,ci ),0 ) << endl;
      lab_partitions( ci, n ) = labels( partitions( n,ci ) ); 
      // getchar();
    }
    cout << endl;
    
    //getchar();*/
  }
  //lab_partitions.print("Final labels: ");
  
}



inline
uvec
k_kmeans::get_id_summ()
{
  uvec indices = sort_index(counts);
  uvec indices_summary;
  
  indices_summary.zeros(N_data);
  
  //partitions.print("Partitions");
  indices.t().print("sorted indices:");
  
  
  double length_summ = floor(N_data*30/100);
  uword j =0;
  uword ne;
  uword pos;
  
  for (uword i = 0; i < counts.n_elem ; ++i)
  {
    pos = indices(i);
    ne = counts( pos ) - 1;
    indices_summary.subvec ( j, j + ne) = partitions(  span(0, ne ) , pos);
    cout << "j= "<< j << " - next j: " << j + ne + 1 << endl;
    j = j + ne + 1;
    uvec tmp = partitions(  span(0, ne ), pos);
    //tmp.t().print("Current indices ");
    //cout << "Current indices " << partitions(  span(0, ne , pos)) << endl;
    //j = j + ne + 1;
    getchar();
    
  }
  
  
  //cout << "p_sum " <<   p_sum << endl;
  indices_summary.resize(length_summ);
  
  //indices_summary.t().print("indices_summary");
  return sort( indices_summary ); // returning only the smaller clusters 
}


inline
uvec
k_kmeans::near_points() // Check it! I implement but I haven't checked
{
  //To select the closest point to the centrioid.
  uvec point_cluster(counts.n_elem);   //Closest point per cluster
  
  vec dist_ci;
  uword  index;
  //double min_val;
  
  for (uword i = 0; i < counts.n_elem ; ++i)
  {
    dist_ci.zeros(counts(i));
    dist_ci = Points_cent( span(0, counts(i) - 1) , i );
    dist_ci.min(index);
    point_cluster (i) = partitions(index,i);
    
  }
  
  
  return point_cluster;
}




inline
umat
k_kmeans::get_partitions()
{
  //partitions.print("partitions");
  return partitions;
  
}






//*****************************************************************************************
//*****************************************************************************************

inline // I think it's more efficient.-- Nov15,2013 
void
k_kmeans::update_kkmeans2()
{
  vec dist(N_cent);
  double dist_cent;
  uword best_ci;
  distortion = 0;
  
  Points_cent.zeros(N_data,N_cent);
  
  new_counts.zeros(N_cent);
  new_partitions.zeros(N_data,N_cent);
  
  
  
  // Precalcular sum2 y almacenar en un vector
  //calcular sum1
  //calcular distance
  
  vec sum2_c;
  sum2_c.zeros(N_cent);
  for (int ce = 0; ce<N_cent; ++ce)
  {
    uvec small_part = partitions( span(0, counts(ce) -1 ), ce );
    sum2_c(ce) = get_sum2(small_part); // in small partition the indices per cluster c are saved
  }
  
  for (int ni = 0; ni < N_data; ++ni)
  {
    for (int c = 0; c < N_cent; ++c)
    {
      double sum1 = 0;

      for (int i = 0; i < counts(c); ++i)
      {
	uword xi = partitions( i , c );
	sum1+=K(ni,xi);
      }  
      
      dist_cent = K(ni,ni) -2*sum1/counts(c) + sum2_c(c) / pow(counts(c),2 );
      dist(c) = dist_cent;
    }
    
    double mind = dist.min(best_ci);
    distortion+=mind;
    //cout << best_ci << ". ";  
    //getchar();
    new_partitions(new_counts(best_ci),best_ci) = ni;
    Points_cent (new_counts(best_ci),best_ci) = mind;
    new_counts(best_ci) += 1;
  }
}

inline
double
k_kmeans::get_sum2(uvec small_part) // -- Nov15,2013 
{
  double sum2 = 0;
  //small_part.t().print("small_part");
  //getchar();
  for (uword i=0; i<small_part.n_elem; ++i)
  {
    uword xi = small_part(i);
    for (uword j=0; j<small_part.n_elem; ++j ) 
    {
      uword xj = small_part(j);
      //cout << K(xi,xj) << " + " << sum2;
      sum2+=K(xi,xj);
      //cout << " = " << sum2 << endl;
      //getchar();
    }
  }
  return sum2;
}